
function [p q] = BDF2(dt)
global U  h  g  p  q   nt dx x h1 p_global ;

dx1=1/dx;
dx2 = 1/(dx*dx);
dt2= 1/dt;

beta =2.0/3.0;
alpha =1.0/3.0;

%---- Construct the spatial Discretization Matrix------%
A = zeros(x,x);
rhs = zeros(2*x,1);
sol=rhs;

A(1,x)=-U*0.5*dx1;
A(1,2) = U*0.5*dx1;
for k=2:x-1
    A(k,k-1)=-U*0.5*dx1;
    A(k,k+1) = U*0.5*dx1;
end


A(x,x-1)=-U*0.5*dx1;
A(x,1)=U*0.5*dx1;

L=zeros(2*x,2*x);
L(1:x,1:x)=A;
L(x+1:2*x,x+1:2*x)=A;

L(1,2*x)   = h*dx2;
L(1,x+2)   = h*dx2;
L(x,x+1)   = h*dx2;
L(x,2*x-1) = h*dx2;

for k=1:x
    L(k,x+k)= -2*h*dx2;
    L(x+k,k)=g;
end

for k=2:x-1
    L(k,x+k-1) = h*dx2;
    L(k,x+k+1)=  h*dx2;
end
I = eye(2*x);  
%------ Build the matrix for the Vertical Structure variable r
b=h*h/3*dx2;

R=zeros(2*x,x);
R(1,x)= b;
R(1,2)= b;
R(1,1)= -2*b;

for k=2:x-1
    R(k,k-1)=b;
    R(k,k+1) =b;
    R (k,k) = -2*b;
end

R(x,x-1)= b;
R(x,1)= b;
R(x,x)= -2*b;

r= zeros(x,1);

% Periodic Boundary Conditions
sol(1:x,1) =p;
sol(x+1:2*x,1) =q;
sol_new =sol;

for n=2:nt+1; 
    % Solve for the vertical structure based on the given values of p and q  
    
    r = R_Solve(sol_new(x+1:2*x,1));
    if (n==2)     
    % Compute RHS which will have contribution from the vertical structure term     
    rhs = (dt2*I - (1-beta)*L)*sol + R*r ;
    sol_new=(dt2*I+beta*L)\rhs;        
    
    elseif (n>2)
       
        rhs = ((dt2*(1+alpha))*I + (alpha - (1-beta))*L)*sol_new - (alpha*dt2)*sol+ (1-alpha)*R*r;   
        sol=sol_new;
        sol_new=(dt2*I+beta*L)\rhs;    
        
    end  

        p = sol_new(1:x,1);
        q= sol_new(x+1:2*x,1);
   time = n*dt;
   
    if rem(time,5)==0
        k=time/5;
        p_global(:,k) = p;
        refreshdata(h1,'caller') % Evaluate p in the function workspace
        drawnow
    end

end
display('Completed Successfully');